Hydrodynamic and hydromagnetic energy spectra from large eddy simulations 
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Direct and large eddy simulations of hydrodynamic and hydromagnetic turbulence have been 
performed in an attempt to isolate artifacts from real and possibly asymptotic features in the 
energy spectra. It is shown that in a hydrodynamic turbulence simulation with a Smagorinsky 
subgrid scale model using 512'' meshpoints two important features of the 4096^ simulation on the 
Earth simulator (Kaneda et al. 2003, Phys. Fluids 15, L21) are reproduced: a correction 
to the inertial range with a k~^^^ Kolmogorov slope and the form of the bottleneck just before 
the dissipative subrange. Furthermore, it is shown that, while a Smagorinsky-type model for the 
induction equation causes an artificial and unacceptable reduction in the dynamo efficiency, hyper- 
resistivity yields good agreement with direct simulations. In the large-scale part of the inertial range, 
an excess of the spectral magnetic energy over the spectral kinetic energy is confirmed. However, a 
trend towards spectral equipartition at smaller scales in the inertial range can be identified. With 
magnetic fields, no explicit bottleneck effect is seen. 

PACS numbers: 44.25.-|-f, 47.27.Eq, 47.27.Gs, 47.27.Qb 



I. INTRODUCTION 

In astrophysical magnetohydrodynamic (MHD) tur- 
bulence, e.g. in stars, accretion discs, the interstellar 
medium, and the intergalactic medium, the magnetic and 
fluid Reynolds numbers are very large. It is therefore 
of great interest to perform simulations with as large a 
Reynolds number as possible. However, the goal of reach- 
ing astrophysical values of the magnetic Reynolds num- 
bers is still far out of reach. The best we can hope for 
is therefore to find asymptotic trends such that one can 
extrapolate into the very large Reynolds number regime. 
However, even that is not really possible as the follow- 
ing estimate shows. As a rule of thumb, for a purely 
hydrodynamical simulation one needs at least an order 
of magnitude for resolving the dissipative subrange, one 
order of magnitude for the bottleneck (a shallower spec- 
trum just before the dissipative subrange), and almost an 
order of magnitude for the forcing to become isotropic. 
This leaves basically nothing for the inertial range-even 
for simulations with 1024^ meshpoints. It is therefore 
only with simulations as big as 4096^ meshpoints [Jjj that 
one begins to see an inertial range. 

In MHD turbulence without imposed field, i.e. when 
the field is self-consistently generated by dynamo action, 
the magnetic energy spectrum peaks at a wavenumber 
that is by a certain factor larger than the wavenumber 
of the kinetic energy spectrum . This factor has been 
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related to the value of the critical magnetic Reynolds 
number for dynamo action, ReM.crit- Specifically, fcmag ~ 

/ckinR.e]^^^j.;^ has been suggested Q, where /cmag and /ckin 
are the wavenumbers of the peaks of the magnetic and 
kinetic energy spectra, respectively, and RcM.crit ~ 35 
This leads to the conclusion that in MHD turbulence one 
needs an even larger Reynolds number than for purely 
hydrodynamical turbulence in order to have a chance to 
see an inertial range. 

What has been found so far is that there is a cer- 
tain range, A:mag ^ A: ^ fed, where the spectral mag- 
netic energy exceeds the spectral kinetic energy 0, 3, 
i.e. there is spectral super-equipartition. While spectral 
super-equipartition is not a priori implausible, it is cu- 
rious that this has not been seen in simulations with an 
imposed field. Such simulations with imposed field have 
recently been performed 0, 0, Q to verify the Goldreich- 
Sridhar theory of MHD turbulence Q . More systematic 
studies of the resulting energy spectra as a function of the 
imposed field strength have been carried out 0, and it 
was found that there is spectral equipartition only when 
the imposed field, Sq, is of equipartition strength, i.e. 
-^0 ^ /^o/Oo^nnsj whcrc /io is the vacuum permeability, 
Po is the mean density, and Urms is the rms velocity. If 
Bq is larger, the magnetic spectrum is always in sub- 
equipartition. 

The case of an imposed field is usually thought to be 
representative of the conditions deep in the inertial range. 
Thus, the observed super-equipartition does seem to be 
in conflict with this result. This is also supported by the 
well known fact that in the solar wind, kinetic and mag- 
netic energy spectra follow a power law with an —5/3 
exponent over several decades |1Q | . In this work we want 
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to elucidate this puzzle by comparing direct simulations 
with simulations using hyperviscosity and hyperresistiv- 
ity, as well as Smagorinsky subgrid scale (SGS) mod- 
elling, in order to imitate larger Reynolds numbers. For 
recent comparisons between direct and Smagorinsky SGS 
simulations; see Refs. [Hill El El, where also decay- 
ing turbulence is considered, albeit only at a resolution 
of 64'^ meshpoints. This was too small to discuss the 
shape of the energy spectra. Recent simulations using 
hyperviscosity have shown that at large enough resolu- 
tion (512'^ meshpoints) the same k~^-^ correction to the 
Kolmogorov k^^^ inertial range spectrum is seen [l^ as 
in the 4096^ meshpoints direct simulations of Kaneda et 
al. . In the present paper we compare these two simu- 
lations also with new Smagorinsky SGS models. 

We need to emphasize that throughout this paper we 
only deal with the case of "non-helical" turbulence, i.e. 
I (w • V X tt) I is negligible (or small compared with ki (it^) , 
where fcf is the typical forcing wavenumber). In some 
sense the case of finite helicity may be regarded as more 
typical (Tij. However, with helicity there is a whole range 
of new problems that need to be addressed. For exam- 
ple, when using hyperresistivity the magnetic field would 
saturate at an artificially enhanced value when there is 
helicity E3 ■ These helicity effects are now fairly well un- 
derstood (see Ref. [3| for a review), but in the present 
paper we discard these complications. 



these expressions are of the general form 



1. 



-Fvisc = -V • (2pj/S) , E,^s = V^J■oJ, (4) 
P 



where 



_ If dui duj 2 J. _ 



(5) 



is the traceless rate of strain tensor. In a direct simulation 
we simply use constant values of j/ and 77, i.e. 



1^ = 1/0, ij = Vo (direct). 



(6) 



In the case of a Smagorinsky SGS model we use i' = I's 
and rj ^ rjs (without constant contributions) where 



5^ '?s = (CmA)VJ2 (Smagorinsky), 

(7) 

where Ck is the Smagorinsky constant. Cm is the mag- 
netic Smagorinsky constant, and A is the filter size, 
which we have set equal to the mesh size. This version 
of the magnetic Smagorinsky SGS model has been stud- 
ied earlier; see, e. g., Ref. nil. Following our experience 
from earlier work |1It| we choose Ck = 0.2, but we vary 
the value of Cm- In simulations with hyperviscosity we 
replace 



II. METHOD 



We solve the compressible non-ideal MHD equations. 



Dm 



1. 



-Vp- 



J X B 



f 



(1) 



where D/Dt — d/dt + it • V is the advective derivative, p 
is the pressure, p is the density, / is an isotropic random 
nonhclical forcing function with power in a narrow band 
of wavenumbers, B is the magnetic field, J = V x B/ fj.Q 
is the current density, and .Fvisc is the viscous force (see 
below). We consider an isothermal gas with constant 
sound speed Cg, so that the pressure is given hy p = c^p 
and p^^'Vp = Cg Vlnp. The density obeys the continuity 
equation. 



Dlnp 
Dt 



= V • u. 



(2) 



The induction equation is solved in terms of the magnetic 
vector potential A, 



dA 

'dt 



U X B Ere 



(3) 



where B — 'V x A is the magnetic flux density, and Ej-cs 
is the electric field due to resistive effects (see below). 

In the following, different combinations of expressions 
for i^visc and £Jics have been explored. In all simulations 



pi/S ^ poi^sV^S, ryJ^r/aVV (hyper), (8) 

in Eq. and use constant coefhcients, referred to as 
v = 1/3 and r] — ri3. Following E3j we use constant 
dynamical hyperviscosity, P0V3 = const, in which case a 
positive viscous heating term can be defined. 

In the present work we only consider cases with 
small Mach number. Compressibility effects are there- 
fore unimportant |l9j . and the continuity equation Q 
can therefore be solved without additional subgrid scale 
terms. We note, however, that by defining suitable aver- 
ages (Favre filtering; see Ref. 0|) the continuity does for- 
mally retain its original form. Likewise, in strongly com- 
pressible flows a turbulent bulk viscosity will be impor- 
tant for smearing out shocks; see, e.g., Ref. [U. Again, 
this is neglected, because we are here only interested in 
nearly incompressible flows. 

It is customary to quote Reynolds numbers based on 
the Taylor microscale A = •\/5wrms/'^ims, where Wrms is 
the rms vorticity, and on the one-dimensional velocity 



dispersion uid, where u^p 



s/3. Hence, we define 



the fluid and magnetic Reynolds numbers for a direct 
numerical simulation as 



Rex 



Re 



(9) 



respectively. Their ratio is the magnetic Prandtl number, 
PrM = v/rj = ReAf/Re, which is unity for all runs. For 
the hyperviscous and Smagorinsky cases we define the 
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Taylor microscale Reynolds number, in analogy to earlier 
work 0, as 

ReA = ReA,o(^)'^', (10) 

where we have defined the effective Kolmogorov 
wavenumber, fcd^eff, whose value is found empirically by 
making the inertial ranges of the spectra overlap as best 
as possible, and Re^.o is a calibration parameter. In an 
earlier paper [l5l | the calibration parameter was found to 
be Rca.o ~ 7.5, which is also the value chosen here. 

We use non-dimensional quantities by measuring 
length in units of 1/ki (where /ci = 27r/L is the small- 
est wavenumber in the box of size L), speed in units of 
the isothermal sound speed Cg , density in units of the ini- 
tially uniform value p = po, and magnetic field in units 

We use periodic boundary conditions in all three direc- 
tions for all variables. This implies that the mass in the 
box is conserved, i.e. (p) = po, where angular brackets 
denote volume averages. We adopt a forcing function / 
of the form 

fix, t) = ^{Nfk(t) exp[ifc(i) • X + i<t>m, (11) 

where x is the position vector, and 5R indicates the 
real part. The wave vector k(t) and the random phase 
—IT < (f){t) < TT change at every time step, so f{x,t) 
is (5-correlated in time. For the time-integrated forcing 
function to be independent of the length of the time step 
St, the normalization factor N has to be proportional 
to 6t^^/^. On dimensional grounds it is chosen to be 
N = /oCs(|A;|cs/i5i)^^^, where /o is a non-dimensional 
forcing amplitude. The value of the coefhcient /o is 
chosen such that the maximum Mach number stays be- 
low about 0.2. Empirically, this is achieved by taking 
/o = 0.02 for all runs discussed below. 

At each timestep we select randomly one of many pos- 
sible wave vectors in a certain range around a given forc- 
ing wavenumber. The average wavenumber is referred to 
as k{. We force the system with nonhelical transversal 
waves, 

/fc = (fe X e) /v/fc2- (fc.e)2, (12) 

where e is an arbitrary unit vector that is real 
and not aligned with k; note that |/fcp = 1. 
For all simulations we use the Pencil Code 
fhttp : //www . nordita . dk/sof tware/pencil-codel 

which is a grid based high order code (sixth order 
in space and third order in time) for solving the 
compressible hydromagnetic equations. 

III. RESULTS 

In an earlier paper [l5j we have shown that hyper- 
viscosity, although it does cause an artificially enhanced 
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FIG. 1; Comparison of energy spectra of the 4096'^ mesh- 
points run l] (solid line) and 512^ meshpoints runs with 
hyperviscosity (dash-dotted line) and Smagorinsky viscos- 
ity (dashed line). (In the hyperviscous simulation we use 
= i-'3 = 5 X 10"^''.) The Taylor microscale Reynolds num- 
ber of the Kaneda simulation is 1201, while the hyperviscous 
simulation of Ref. [lfll | has an approximate Taylor microscale 
Reynolds number of 340 < ReA < 730. For the Smagorinsky 
simulation the value of ReA is slightly smaller. 

bottleneck effect in purely hydrodynamic turbulence, it 
does not affect the inertial range if the resolution is 
large enough. Instead, hyperviscous simulations with 
512'^ meshpoints reproduce the k^^-^ correction with 
wavenumber k. This was first found by Kaneda et al. 
Il|. We begin by comparing these results with simula- 
tions where Smagorinsky SGS viscosity is used. 



A. Hydrodynamic turbulence 

In Fig. ^ we compare kinetic energy spectra of runs us- 
ing ordinary viscosity (4096'^ meshpoints, solid line) by 
Kaneda et al. with runs using Smagorinsky viscosity 
(512'^ meshpoints, dashed line) and runs using hyper- 
viscosity (512^ meshpoints, dash-dotted line). Since the 
simulation with 4096'^ meshpoints and ordinary viscosity 
is the largest direct simulation to date, we use it as our 
benchmark. The spectra for the runs with hyperviscos- 
ity and Smagorinsky viscosity have been scaled by em- 
pirically determined factors 1.1 and 0.95, respectively, so 
as to make the spectra overlap within the inertial range. 
However, these scaling factors are still well within the 
range over which the spectra fluctuate in time. 

We see that at all scales (including those of the bot- 
tleneck) the simulation with Smagorinsky SGS modeling 
is surprisingly similar to the benchmark result. Further- 
more we see that at large scales and in the inertial range 
the run with hyperviscosity agrees well with the bench- 
mark result. The bottleneck is however greatly exagger- 
ated in height, even though the width is the same jl5|. 
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Most important is perhaps the k~°-^ correction to the 
usual k"^/^ inertial range scaUng. The same correction 
is seen in ah three simulations. The k~^'^ correction 
is stronger than the usual intermittency correction pre- 
dicted by the She-Leveque model J2J, which would only 
predict a fc~0 03 correction. This strong correction may 
be an artifact of the absence of a well resolved subinertial 
range j2^]. This would be in some ways just opposite to 
the emergence of a shallower spectrum near the dissipa- 
tive cutoff wavenumber if the dissipative subrange is not 
well resolved p3 |. 

The only major discrepancy between the Smagorinsky 
and direct simulations is the lack of a sharp decline of 
the spectral energy toward the right of the bottleneck. 
In order to understand this difference, we must first of 
all recall that our Smagorinsky simulation did not have 
any explicit (constant) component at all (z/q = 0). There- 
fore, if the Smagorinsky model was a perfect subgrid scale 
model, it should represent the infinite Reynolds number 
case. The bottleneck would then be far to the right and 
outside the graph, so one should only have a pure Kol- 
mogorov spectrum. The reason for the bottleneck in the 
Smagorinsky case is therefore related to the fact that we 
are still working here with an ordinary diffusion operator 
using just a variable viscosity coefficient. Therefore, the 
standard explanation for the bottleneck still applies; it is 
caused by strongly nonlocal interactions in wavenumber 
space, corresponding to wave vectors forming strongly 
elongated triangles. Close to the viscous cutoff wavenum- 
ber, these interactions prevent the disposal of energy 
from the end of the inertial range, which then causes the 
pileup of energy near the dissipation wavenumber |24| . 
The same argument also applies to the current case of 
Smagorinsky viscosity. In conclusion, the reason for the 
discrepancy between direct and Smagorinsky simulations 
to the right of the bottleneck is that the Smagorinsky 
model tries to maintain pure Kolmogorov scaling every- 
where, but fails to do so just before the cutoff wavenum- 
ber imposed by the finite mesh resolution. 



TABLE I: Summary of the most important runs. The mean- 
ing of entries in the columns for u and 77 depends on the entry 
for 'Method', as explained in the text. In the Smagorinsky 
cases ordinary viscosity is neglected, i.e. u — 0. Except for 
Method O, the resulting values of kd^cB, and hence also of 
Rca, are uncertain within ~ 40% 



Run 


Res. 


Method 


V 




■n 






ReA 


A 


1024'' 





8 X 10- 


-5 


8 X 10" 


- b 


143 


200" 


Bl 


128^ 


II 







1 X 10" 


-9 


180 


180 


B2 


256=^ 


II 







3 X 10" 


11 


330 


270 


B3 


512^ 


II 







5 X 10- 


13 


700 


450 


CI 


128^ 


III 


1 X 10" 


-9 


1 X 10" 


-9 


180 


180 


C2 


256^ 


III 


3 X 10" 


11 


3 X 10" 


11 


330 


270 


C3 


512^ 


III 


5 X 10" 


13 


5 X 10" 


13 


700 


450 



"Note that in Ref. 25] the value of ReA was based on the 3- 
dimensional velocity dispersion, so the non-magnetic equivalent of 
Run A was quoted with Rca = 350. 



Smagorinsky, 128^ 
Smagorinsky, 64^ 
Ordinary, 10243 ^^.^ 
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FIG. 2: Total magnetic and kinetic energies for runs 
with 128'^ (solid line) and 64'^ (dashed line) meshpoints and 
Smagorinsky diffusion and resistivity (Method I) compared 
with a direct simulation with 1024^ meshpoints (Method O, 
horizontal dotted lines) . Note the lack of convergence for any 
value of Cm. 




B. Hydromagnetic turbulence 

For the MHD case we use a 1024'^ meshpoints sim- 
ulation with ordinary viscosity as our benchmark 
We compare with the SGS model where Smagorinsky 
schemes are used both for the velocity and the magnetic 
fields. In the following we refer to this as Method I. 
We also compare with cases where we use hyperresis- 
tivity. In the momentum equation we use either the 
usual Smagorinsky SGS model, which is referred to as 
Method II, or we use hyperviscosity (Method III). The 
results of these three methods are compared with those of 
direct simulations (Method O). In summary, the different 



methods considered here are 

Method I: va, and rys (full Smagorinsky), 

Method II: v<g, and 773 (Smagorinsky/hyper), 

Method III: and 773 (full hyper). 

Method O: vq and 770 (benchmark). 

We have listed the relevant runs in Table 

In Fig. 121 we show that the agreement between the 
results of Method I and the benchmark is poor. The 
dynamo-generated magnetic energy remains far below 
the benchmark target. The largest value of the mag- 
netic energy is reached for Cm = 0.3, but even then it is 
only about 30% of the target value. 

In order to understand the reason for the poor perfor- 
mance of the magnetic Smagorinsky model (Method I) 
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FIG. 3: Comparison of magnetic and kinetic energy spectra 
of runs using Method I with 128'^ meshpoints and various 
values of Cm- 



we plot in Fig. |31 kinetic and magnetic energy spectra for 
various values of Cm- Clearly, for Cm < 0.3 both kinetic 
and magnetic spectra diverge toward large wavenumbers. 
This shows that this model becomes unphysical and can- 
not be used for too small values of Cm- For Cm = 0.5, 
on the other hand, magnetic and kinetic spectra fall off 
at large wavenumbers. However, the effective resistivity 
of the magnetic Smagorinsky scheme is apparently too 
large for Cm = 0.5, so that the dynamo is suppressed. 
The poor performance of this model is not too surprising 
if one recalls that it is a rather crude method in that it 
deals with the small scales only in a diffusive manner. 
We also note that the Smagorinsky SGS model has, to 
our knowledge, never before been tested in the context of 
dynamo action. We conclude that using the Smagorinsky 
SGS model for the magnetic field docs not give satisfac- 
tory results. Therefore, from now on, we discard it as 
inappropriate for our purpose. 

We see from Fig.0]that the compensated spectra with 
only 128'^ meshpoints, using Methods II and III, match 
the benchmark quite well at all scales down to the dissi- 
pative scale. We have compensated the energy spectra by 
/c^/'^erp^^'^, such that a Kolmogorov-like spectrum would 
appear flat. Here ex = ek + fM, where ck and cm are 
the kinetic and magnetic dissipation rates, respectively. 
The kinetic energy spectrum of the 1024"^ run has how- 
ever been multiplied by 1.3 in order to make all spectra 
overlap. We believe the shift is due to the fact that the 
1024'^ run has not been run for very long, and the aver- 
age dissipation rate, ex, has not yet fully converged, even 
though the slope converges generally much quicker. From 
the general agreement between the three runs shown in 
Fig.^we conclude that, for our purpose, Methods II and 
III give useful results. 

In Fig. [51 we compare compensated spectra for three 
simulations which all use Smagorinsky viscosity and hy- 



FIG. 4: Comparison of magnetic and kinetic energy spectra 
of runs with 1024'^ meshpoints and normal diffusion (Run A, 
solid line) with 128'^ meshpoints and hyperdiffusion (Run CI, 
dash-dotted line), and with 128^ meshpoints and Smagorinsky 
viscosity and hyperresistivity (Run Bl, dashed line). Note 
that both the magnetic and kinetic energy spectra for the 
three runs are very similar for fc/fcd < 0.1. 




0.01 0.10 

^/^d.eff 



FIG. 5: Magnetic and kinetic energy spectra for runs with 
128^ (Run Bl), 256^ (Run B2) and 512^ (Run B3) meshpoints 
where all of them use Smagorinsky viscosity and hyperresis- 
tivity (Method II). Note the approach of the kinetic energy 
spectra towards the magnetic energy spectra at a point that 
is well before entering the bottleneck and the dissipative sub- 
range. 



perresistivity, but have different Reynolds numbers. We 
see that, unlike the purely hydrodynamic case, the dis- 
sipative subranges do not collapse onto the same func- 
tional form for different Reynolds numbers. On the other 
hand, for purely hydrodynamical simulations |15| the dis- 
sipative subranges collapse very well onto the same func- 
tional form and the inertial range simply becomes longer 
for larger Reynolds numbers. Furthermore, in Fig 1 of 
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FIG. 6: Magnetic and kinetic energy spectra for runs with 
128^ (Run CI), 256^ (Run C2) and 512^ (Run C3) mesh- 
points where all of them use hyperviscosity and hyperresis- 
tivity (Method III). 




0.01 I J 

0.001 0.010 0.100 1.000 

^/^d.eff 

FIG. 8: Sketch of kinetic and magnetic energy spectra, follow- 
ing the Miiller and Grappin phenomenology. Note the slight 
super-equipartition just to the right of the peak of Euik) and 
the asymptotic equipartition for large wavenumbers. 



1.00 



0.01 







'/ 

■ A 
1 


w * 
w ^ 

Hyperdiffusion \\ 
_ Smagorinsky \ 
V 



0.01 0.10 

^/^d.eff 



FIG. 7: Magnetic and kinetic energy spectra for runs with 
512^ meshpoints and hyperviscosity and hyperresistivity (Run 
C3, solid line) and Smagorinsky viscosity and hyperresistivity 
(Run B3, dashed line). Note the mutual approach of kinetic 
and magnetic energy spectra before entering the dissipative 
subrange. 

Ref. ^3 we see that the bottleneck is similar and con- 
stant for all Reynolds numbers. Again, in the MHD sim- 
ulation we see nothing similar. 

In Fig. 13 we have shown the same as in Fig. |S1 but us- 
ing hyperviscosity instead of Smagorinsky viscosity. We 
clearly see that the tendency is the same in both fig- 
ures. Since the bottleneck effect is quite different for pure 
hydrodynamical simulations with Smagorinsky viscosity 
and with hyperviscosity (see Fig. QJ, it is reasonable to 
assume that the tendency we see is robust and not due 
to the specific modeling applied, but that it is a physical 
effect. 



Finally, we compare in Fig. [T] spectra of Smagorin- 
sky and hyperviscous simulations using the highest avail- 
able resolution of 512'^ meshpoints. Again, note that 
the spectra for hyperviscous simulations and those with 
Smagorinsky SGS modeling are almost identical. Fur- 
thermore, there is no range where both kinetic and mag- 
netic energy spectra are parallel. Together with the re- 
sults of Figures [3 and El we therefore conclude that we 
have not yet reached Reynolds numbers large enough to 
show an inertial range. 

IV. SPECULATIONS ON ASYMPTOTICS 

The direct MHD simulations of Ref. 2] have sug- 
gested the presence of a super-equipartition range where 
Euik) ~ 2.5i?K(fc)- However, the spectra still showed 
some weak bending, indicating that a proper inertial 
range has not been reached even at a resolution of 1024^ 
meshpoints ,26] . The present SGS models reproduce 
the spectral super-equipartition of magnetic over kinetic 
spectral energy (Fig. [7|), but they also show now more 
clearly that the two spectra are not parallel to each other. 
Instead, they approach each other in such a way that the 
compensated kinetic energy spectrum shows a strong up- 
rise. 

One might argue that the uprise at the end of the com- 
pensated kinetic spectrum is just a strong bottleneck. 
This is however unlikely since both SGS models give the 
same uprise, even though in purely hydrodynamic tur- 
bulence the hyperviscosity model is known to produce 
a much higher bottleneck than the Smagorinsky model 
(Sec. IhT^ . Furthermore, in hydrodynamic turbulence 
the width of the bottleneck is independent of Reynolds 
number, whereas in the present case it appears to become 
wider with increasing Reynolds number. This suggests 
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FIG. 9: Residual spectrum, En ~ Em — Ek, compensated 
by for the same runs shown in Fig. Q Negative values 

of -Em — Ek are indicated by dotted and dash-dotted lines for 
hyperdiffusion and Smagorinsky runs, respectively. 

that the uprise in the MHD case is a true large scale fea- 
ture of the spectrum, and independent of the dissipative 
subrange. 

Next, we recall that in simulations with an imposed 
magnetic field, the magnetic and kinetic energy spectra 
are found to be in approximate equipartition only when 
the field strength is of the order of B^q . Such simula- 
tions are thought to be representative of the small scale 
end of the inertial range of any MHD simulation, even if 
the field is generated by a small scale dynamo as in the 
present case. Assuming that this interpretation is cor- 
rect, it would support our previous suggestion that the 
spectral super-equipartition was only a non-asymptotic 
feature, confined to the large scales, and not a true in- 
ertial range feature. We are therefore led to believe that 
for much larger Reynolds numbers the kinetic and mag- 
netic energy spectra might converge. Qualitatively, this 
can be reproduced by the phenomenology proposed re- 
cently by Miiller and Grappin [23, 13 ■ According to 
their theory, the total energy Et = Em + Ek follows still 
the expected k~^/^ spectrum, while the residual energy 
E^ = \Em — Ek\ follows a k~'^/^ spectrum. In Fig.|Slwe 
produce such an example with 

for fc > 1 (in arbitrary units). Using the fact that in the 
inertial range Em exceeds Ek by about a factor a — 2, 
we reconstruct Em and introduce an additional k^ subin- 
ertial range, so we write 

EM{k) - ^[Erik] + Enik)]/[1 + {k/kM)-^^/\ (14) 

with fcM — 5. The kinetic energy obtained by assuming 
that the total energy is constant is 

EK{k) = [Erik) - EM{k)]/[l + (fc/fcx)""/'], (15) 



where we have included a different subinertial range be- 
low kK — 1.5. The resulting spectra shown in Fig. |H| 
reproduce surprisingly well the basic features suggested 
by our SGS simulations of Fig. [7| 

In order to see how well our simulations reproduce the 
anticipated scaling of the residual spectrum, we 

plot in Fig. |5| the appropriately compensated E'r spec- 
trum. Clearly, the residual spectrum is still curved, but 
it remains reasonably straight within about half an order 
of magnitude in wavenumbers. 

V. CONCLUSIONS 

The results of subgrid scale models should always be 
taken with great care. Even if their results can be trusted 
in one case (e.g. in the case without magnetic fields), they 
may not give reliable results in another case (e.g. in the 
presence of magnetic fields and dynamo action). How- 
ever, once we begin to see detailed agreement between 
SGS models and direct simulations, it may be possible 
to use this agreement to justify the use of the SGS model 
in more extreme parameter regimes that are currently 
inaccessible to direct simulations. 

In the present work we have shown that the Smagorin- 
sky SGS model with a resolution of 512'^ meshpoints is 
able to reproduce the hydrodynamic turbulence spectra 
of a direct simulation at an almost 10 times larger reso- 
lution (Fig. |21l . On the other hand, an extension of this 
model to the MHD case with dynamo action leads to ob- 
vious problems (the intensity of the dynamo is artificially 
suppressed). However, using hyperresistivity instead of 
a Smagorinsky-type SGS model leads to fair agreement 
between the 128'^ SGS simulation and the nearly 10 times 
larger direct simulation (Fig.^. Thus, having validated 
the SGS model at 128^ meshpoints, we may be justified 
in proceeding further to a resolution of 512^ meshpoints 
(Fig. 01- Here, a new and yet unconfirmed feature arises: 
a tendency towards spectral equipartition. This, together 
with the knowledge that there is spectral equipartition 
with imposed fields of equipartition strength |9| , suggests 
a spectrum that might look like what is shown in Fig. |S1 

Obviously, we will not be able to verify this result in 
the immediate future. Although it may soon be possible 
to obtain the resources necessary to do a 4096'^ MHD sim- 
ulation to validate the results of Fig. [7| yet another order 
of magnitude in improved resolution will be necessary to 
test the hypothesis sketched in Fig. |S1 Our results may 
therefore serve as a justification for using future comput- 
ing resources for this type of problem. 
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